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Abstract 

We present a computer-assisted approach to locating approximate coarse optimal switching policies be- 
tween stationary states of chemically reacting systems described by microscopic/stochastic evolution rules. 
The "coarse time-stepper" constitutes a bridge between the underlying kinetic Monte Carlo simulation and 
traditional, continuum numerical optimization techniques formulated in discrete time. The approach is il- 
lustrated through two simple catalytic surface reaction models, implemented through kinetic Monte Carlo: 
NO reduction on Pt, and CO oxidation on Pt. The objective sought in both cases is to switch between two 
coexisting stable stationary states by minimal manipulation of a macroscopic system parameter. 

keywords: equation-free, dynamic optimization, coarse timesteppers, kinetic Monte-Carlo, Hooke- Jeeves 

1 Introduction 

The computation of operating conditions that are optimal with respect to a predetermined objective (derived 
from economic considerations, safety requirements and product quality constraints) has been, for many decades, 
an exciting research area for chemically reacting systems. Increase in computational power has led to increased 
attention to dynamically evolving processes and the search for optimal time- varying operation protocols, both for 
lumped and spatially distributed systems (e.g., [47, 49, 1, 46, 3]). Interestingly, as sensing and actuation become 
increasingly more resolved in space and time for emerging chemical processes, spatiotemporally complicated 
operating policies can be considered (e.g., [50, 37, 19]). The problem of computing optimal transition paths over 
complicated free energy surfaces (see e.g., [9]) is also computationally related, especially considering Jarzynski's 
relation between equilibrium free energy and nonequilibrium work (see e.g., [44]). 

Existing computational approaches for solving dynamic optimization problems at the deterministic, contin- 
uum level may involve (a) formulation as a temporally discretized problem (both for the process and for the 
operating variable(s)) simultaneously, and solution using large sparse linear algebra techniques (e.g [48, 5, 47]); 
(b) formulations involving direct integration of the model equations in time, keeping track of possible constraint 
violations and temporal discretization of the operating variables [8, 10, 49, 2]; or possibly (c) using efficient so- 
lution algorithms within a dynamic programming formulation [29, 31, 4]. Knowledge of a macroscopic process 
model, in the form of macroscopic mass balances closed through appropriate constitutive expressions -such as 
chemical kinetic rate formulas-, is a fundamental prerequisite for these computational solution strategies. 

In contemporary engineering applications, however, we are often faced with problems for which the available 
physical description is in the form of atomistic / stochastic evolution rules (e.g., kinetic Monte Carlo, Lattice 
Boltzmann, molecular dynamics, Brownian dynamics) while the design, optimization or control is required at 
a coarse-grained, macroscopic level. The lack of a closed-form process model can in principle be circumvented 
during the formulation of steady state optimization problems through the construction of identification algorithms 
"wrapped" around black box simulators (e.g., see [41]). Through design of computational experiments, sampling 
and estimation, one can determine the form of "surrogate" models that are subsequently used to identify local 
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or global (e.g., see [22, 35]) optima. In another approach, near optimal solutions are found for unconstrained 
optimization problems, under the assumption that absolute and relative error bounds are known for the computed 
objective function values [24]. 

Over the last few years we have been developing an equation-free computational approach enabling microscopic 
/ stochastic simulators to directly perform system-level tasks, such as coarse integration, stability analysis, bi- 
furcation/continuation and feedback control, thus circumventing the derivation of explicit macroscopic evolution 
equations [45, 14, 25, 40]. In this work we demonstrate the applicability of this computational enabling technology 
to coarse-grained optimization tasks; in particular, we present a formulation methodology capable of addressing 
dynamically evolving processes, to derive optimization formulations that can be solved employing standard, off- 
the-shelf direct search algorithms. The approach is applied to computationally identify coarse optimal operating 
parameter policies capable of switching the (expected) behavior of kinetic Monte Carlo simulators from one sta- 
tionary state to another (alternatively, from the bottom of one potential well to the bottom of another); this is 
illustrated through kinetic Monte Carlo (KMC) simulations of two simplified models of heterogeneous catalytic 
chemical reactions. Both these systems are characterized by two "coarse-grained" stable stationary states. We 
seek optimal (for a particular definition of a cost function) parameter variation policies that will switch the KMC 
simulation from one stationary state to the other within a finite time interval. 

The paper is organized as follows: In the following Section we present the problem we want to solve and our 
two illustrative examples. We then present our formulation for coarse-grained computational optimization. In 
the following two Sections we present detailed computational results for each of our illustrative problems; we then 
conclude with a brief discussion. 



2 Process description 

We investigate microscopic/stochastic processes for which we believe that the coarse-grained, expected dynamics 
can be well approximated in closed form by an equation of the general type 

x = F{x,p) (1) 

but where the right-hand-side of the evolution law, F, is not available in closed form. Here, x{t) e IR" is a state 

dx 

variable vector, t is the time, x is the time derivative of a;, — , and p e "P™ is the vector of process parameters 

at 

(j)m (- jj^rn jg ^]^g subset of acccssibfe values of the process parameters). The state variables x{t) in the model 
of Eq.l are coarse-grained observables of the microscopic simulation - typically a few lower order moments of 
an atomistically or stochastically evolving distribution (e.g., a concentration, or, in our heterogeneous catalytic 
reaction examples, a surface coverage, the zeroth moment of the distribution of adsorbates on the surface). The 
(unavailable) equation for the expected behavior of the process may possess, at fixed process parameter values, 
one or more steady states Xss,i- We will assume that initially the microscopic process is in the neighborhood of 
one of these stable, coarse-grained stationary states, represented by Xss,i for the system of Eq.l. 

2.1 The coarse time-stepper 

The basis of our approach is the computation of a deterministic optimal policy -a time-dependent protocol for 
changing one of the operating parameters of the process- so as to attain a particular goal for the expected 
dynamics of the process (the "coarse-grained dynamics"): switching from one stationary state to another one. 
We want to accomplish this by acting directly on the microscopic simulator, thus circumventing the necessity 
of first deriving a closed form evolution equation for these coarse-grained dynamics. This deterministic policy 
for the coarse-grained behavior will then be applied to individual realizations of the process. As we will discuss 
below, it is convenient in our approach to reformulate the coarse dynamics in discrete rather than continuous 
time. The coarse-grained evolution law then takes the form: 

Xi+i = GT{xi,Pi+i) (2) 

where ti is the time at the beginning of i-th time interval, T is the time interval duration (reporting horizon), Gt 
represents the evolution of Eq.l for constant process parameter value p{t) = Pi+i, € {ti,ti + T], initialized at 
Xi and evolved for time T, arriving at state a;,+i. 
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Conventional algorithms for the solution of optimization problems involving Gt incorporate frequent calls 
to a subroutine that evaluates Gt and/or the action of its derivatives in phase and parameter space on initial 
conditions. Equation- free based algorithms estimate the same quantities by short bursts of (possibly ensembles of) 
microscopic simulations conditioned on the same macroscopic initial conditions ([45, 14, 25, 40, 26]). The coarse 
time-stepper constitutes such an estimate of the discrete-time, macroscopic input-output map Gt obtained via 
the kinetic Monte Carlo simulator. Through a lifting operator the macroscopic initial condition is translated into 
several consistent microscopic initial conditions (distributions conditioned on a few of their lower moments). In 
some cases constructing microscopic initial conditions consistent with macroscopic observable quantities is easy 
(e.g., constructing distributions with prescribed means and variances); in other cases (e.g., when wc want to 
prescribe pair probabilities on a lattice, or even a pair correlation function) an optimization problem may need 
to be solved to successfully lift; short runs with constrained dynamics algorithms ([39, 13, 12]) can also help in 
preparing such consistent microscopic initial conditions. 

This ensemble of microscopic initial conditions is then evolved microscopically, in an easily parallelizable 
fashion (one consistent realization per CPU) for a relatively short time. The results are averaged through a 
restriction operator back to a macroscopic "output" ; it is precisely this output that traditional algorithms simply 
compute through function evaluations when the evolution equations are available in closed form. As extensively 
discussed in [33, 25] part of the microscopic evolution is spent in a "healing" process - the higher moments, 
which have been initialized "wrong" quickly relax to functionals of the low order moments (our state variables). 
A separation of time scales (fast relaxation of the high moments to functionals of the low ones, and slow - 
deterministic- evolution of the low ones) underpins the existence of deterministic coarse grained evolution laws. 
The dynamics of the evolving microscopic distribution moments in the problems we study constitute thus a 
singularly perturbed system. The requirement of finite time microscopic evolution (necessary for the moment 
healing process) conforms with the discrete-time formulation of the coarse optimization problem, which is common 
in many optimization algorithms (see section 3 below). 



2.2 Numerical experiments 

Wc illustrate the proposed combination of traditional optimization techniques with stochastic simulators through 
KMC realizations (using the stochastic simulation algorithm, proposed by Gillespie [15, 16, 17]) of kinetic models 
describing catalytic NO reduction by H2 and CO oxidation by O2, respectively. 

We initially apply the proposed method to a drastically simplified kinetic model of NO reduction by H2 on 
Pt surfaces; the model involves Langmuir adsorption, first order dcsorption, and chemical reaction requiring two 
neighboring vacant sites. The mean field Langmuir-Hinshelwood approximation of the kinetic mechanism we will 
model consists of a single Ordinary Differential Equation for the surface coverage of NO: 

ja 

^ = a{i-e)--ie-k{i-efe. (3) 

Here 6 describes the surface coverage of adsorbed NO, a, 7 are the NO adsorption rate (incorporating the 
dependence of gas phase pressure of NO) and desorption rate constants respectively, and k is the reaction rate 
constant. The reaction term is third order due to the need for two free adjacent sites for the adsorption of H2- In 
Figure 2 we present the deterministic bifurcation diagram in the form of coverage Ogg at steady-state as a function 
of fc for a = 1 and 7 = 0.01. We observe a range of k values for which the system exhibits multiple steady-states; 
the higher and lower ones are locally stable, while the middle one is unstable. 

The proposed approach is also applied to a kinetic Monte Carlo description of the CO oxidation reaction 
mechanism, using a simplified model for the reaction kinetics of the form A -\- I/2B2 — > AB. The kinetics here 
involve Langmuir adsorption for A, dissociative adsorption for B2, and a second order surface reaction whose 
products are immediately desorbed. The mean field approximation equations for this mechanism in the absence 
of adsorbate interactions would consist of a set of two ODEs [33]: 

^ = a(l - - Ob) - 7^a - ^KOaOb 

(4) 

^ = 2/3(1 -Oa- - AkrOAOB 

where 6a and Ob describe the surface coverage of adsorbed CO and O2 respectively, a and (} denote the adsorption 
rate constants of CO and O2 respectively, 7 the desorption rate constant of CO, and with kr we denote the 
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oxidation rate constant. In Figure 3 we present the deterministic bifurcation diagram representing the coverage 
of A, 9a, and B, 6b, respectively, as a function of the adsorption rate constant /3 for values of a = 1.6, 7 = 0.04 
and kr — 1. Wc observe that the system exhibits multiple steady states, the higher and lower ones being locally 
stable and the middle one unstable. Note that when the steady-state value of A coverage is low, then the coverage 
of B is high and vice- versa. A description of the use of the coarse KMC time-stepper in obtaining "coarse" versions 
of this bifurcation diagram may be found in [33]. 



3 Coarse computational optimization 

Computing the temporal profiles of the process parameters that cause the transition of the coarse-grained system 
from one initial stationary state to a different final one can be formulated as an optimization problem; the 
objective is to minimize an integral cost function over time: 



mm 



/ Q{t,x,p)dt (5) 
Vo 



where Q is a continuous scalar cost function. The constraints for this coarse optimization problem are the 
(unavailable) coarse process evolution equations Eq.l, the initialization at one of the coarse stationary states for 
P = Pss, the requirement of termination at a different stationary state for the same process parameter value Pss, 
as well as, possibly, other inequality constraints g{x,p) (e.g., that surface coverages should remain positive and 
sum up to less than 1 at all times): 

x-F{x{t),p{t))=0, g{x,p)<0 
p{0)=Pss, liinp{t)=pss (6) 

a;(0) = Xss,i, lim x{t) = Xss.2- 

t — *oo 

This constitutes an infinite dimensional problem in continuous time. Direct solution methods are based on the 
calculus of variations. Semi- infinite programming approaches provide us with the necessary mathematical tools 
to solve such problems with finite time horizon [20] , through discretization of the temporal domain. 

Another approach consists of approximating this problem through a finite time horizon problem with a final 
state penalty, which is (in our case) solved in discrete time. This results in a finite dimensional, generally nonlinear, 
optimization problem which -if the coarse equations are available- could be solved using available optimization 
techniques. For example, discretizing the process time [0,tf] in A'' time intervals of length T (not necessarily 
constant) and assuming that the process parameters remain constant within each interval, results in the following 
optimization program with (A?' -|- 1) x (n -|- m) variables and nx {N + 1) +m + n equality constraints: 

N 



min / Qd{t,Xi,Xi_i,Pi)dt 



+W{R{\x{tf)-x,s,2\-e)) 

s.t. (7) 
Xi = GT{xi-i,pi) = Q, i = l,2,...,N 

9d{X0,---,XN,P0,---,PN) < 0, 
PO=Psa, Xq = Xss,1- 

Here Qd is analogous to the Q function for the discrete values of the state, gd is similarly analogous to g, W(-) 
is a class K scalar function, R{-) is the ramp function and e is a value for which |xAr — 2I < e <=> W = 0. 
Appropriate final state penalty contributions to the objective function take the place of final state constraints 
at infinite time as stated in the original formulation of the problem; the final state is restricted to be (in finite 
time) within a neighborhood of the final stationary state. This fully discrete-time formulation is ideally suited for 
linking with a coarse time-stepper. It is also in principle possible to discretize only the process parameter temporal 
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behavior, resulting in the following formulation (coined control vector parameterization [48]) with N x m + 2n 
variables and n equality constraints: 

min / Q{t,x,p)dt 

+W{R{\x{tf)-XssM-e)) 

s.t. (8) 
a; - ^'(a;,^) = 0, g{x,p)<0 

^ t 1 

Pit) = Ep«n(- - ^ + 2) +PssHi-t) 

i=l 

where n(-) is the standard boxcar function and H{-) denotes the Hcavisidc function. In cases where the explicit 
form of Eq.l is unavailable, the state evolution is provided through direct simulation of the system. 



3.1 Solution methodology 

Traditional discrete time optimization schemes would repeatedly call, during the solution process, a numerical 

integration subroutine for the system equations (and possibly variational and sensitivity integrations for the 
estimation of derivatives with respect to state variables and parameters). This call is now substituted by the 
coarse time-stepper; the most important numerical issue is that of noise, inherent in the lifting process and 
the stochastic simulations, and the variance reduction necessary to estimate the state or its various derivatives. 
Simulations of different physical size systems (different lattice sizes in our simulation) are characterized by different 
values of variance. For the type of simulations in this paper, changing the physical size of the simulated domain 
on the one hand, and changing the number of copies of the simulation on the other, have comparable effects in 
reducing the variance of the simulation output; this, however, is not generally the case in KMC simulations. When 
a switching policy for a particular physical size system is required (e.g., for nanoscopic reacting systems, such as 
the chemical oscillations on Field Emitter tips [43]) variance reduction can be affected through a larger ensemble 
of consistent microscopic initializations (see also [34, 36] for variance reduction techniques). Simulation noise 
affects both function evaluations and (numerical) coarse derivative evaluations, and thus becomes an important 
element of the approach. It is interesting, however, to observe that massively parallel computation can reduce 
the wall clock time required for the computation distributing different microscopic initializations/realizations to 
different CPUs. 

Due to the presence of noise in the coarse time-stepper results, the use of optimization algorithms that are 
specifically designed to be insensitive to noise becomes desirable and even necessary (the numerical estimation of 
derivatives is, of course, highly susceptible to noise). A class of search algorithms that fulfill this requirement are 
ones that use only function evaluations to search for the optimum, such as the Luus- Jaakola [32] , Nelder-Mcad and 
Hooke-Jeeves algorithms. Algorithms that use local and bounded approximations of the Jacobian matrix of the 
objective function with respect to the process parameters have also been developed, such as the implicit filtering 
algorithm; the reader may refer to [23, 27] for a review of these methods. In our approach, an optimization 
"wrapper" is constructed around an on demand estimate of the system coarse-grained model, constructed by 
computational experimentation with the microscopic solver. There is a close relation to the algorithms in [23, 
27, 24], where comparable wrappers are constructed around continuum models that are difficult or expensive 
to evaluate (e.g., the ones that come from solutions of large PDEs). It is also remarkable that this "wrapper" 
approach, which performs the equation-free solution of an optimization problem through appropriate design of 
computational experiments, can be also in principle be used for the equation-free optimization of experimental 
systems. Indeed, if enough control authority exists for experiments to be initialized in detail, these optimization 
"wrappers" become protocols for the design of laboratory experiments leading to an optimum [18, 30, 38]. 

In all the above iterative algorithms [23], an appropriate initial guess of the process policy profile and a 
set of search directions for the m x N coarse-grained variables is required. The usual set of search directions 
(also used in our numerical experiment) is the unitary basis for JR"^^^ . Specifically for stencil-based search 
algorithms, important parameters also include a vector, the elements of which are the maximum distances the 
algorithm should venture from the current position during the new direction search step at each iteration. These 
perturbation distances are called scales. 



5 



Selection of these scales in the search algorithm requires estimates of noise bounds. A scale that is "too 
small" not only leads to increased computation time with no apparent advantages, but, especially in the case of 

algorithms that compute approximations of the Jacobian, can lead to grossly erroneous results for the search. 
The following time-stepper protocol yields, in our case, simulation results of bounded ensemble variance: 

1. Initialize time-stepper. Set: 

• lattice size Ni, 

• solution replica ensemble size Mr, 

• minimum Mmin and maximum M^aax- 

• ensemble variance limit dmax and 

2. In each time step, the time-stepper: 

(a) simulates system for the desired parameter value Mr times (this can be affected through different 
microscopic initializations and/or different random seeds in the Monte Carlo process). 

(b) Computes the standard deviation divided by the average of the sample, d. 

(c) If > dmax and Mr < Mmax, adjusts Mr > Mmin, repeat step (a). 

Adaptively adjusting the number of realization in the ensemble can thus help in the choice of appropriate scales. 
The lattice size Ni of the time-stepper routine affects the expected evolution of the process observables; for the 
SSA simulations in this paper, as the size of the simulation increases the expected behavior of the stochastic 
microscopic process becomes closer and closer approximated by the mean-field description represented by the 
ODE system [16]. In our computations in this paper, the lattice size of our time-stepper is chosen so that 
the discrepancy between the expected behavior obtained from the KMC simulations and the ODE solutions is 
negligible [11]; this way we can validate the directly computed optimal switching policies by comparing them 
to the mean-field, deterministic ones. Changing the lattice size will, in general, affect the expected behavior of 
a stochastic process; the variance of the time-stepper results will also be affected. We will use multiple replica 
simulation runs M,. (with the same process settings and macroscopic initial conditions) to control fluctuations in 
the estimation of the expected results. 

3.2 Direct search method: The Hooke-Jeeves algorithm 

The processes that we investigate often lead to optimization problems, where the objective functions exhibit a 
large number of small local optima that in effect "hide" the directions along which the objective function decreases. 
This is due to the stochastic nature of the microscale simulations that are used to infer the dynamic behavior of 
the coarse variables. 

A class of search algorithms that avails itself to the solution of such problems is direct search algorithms. 
These methods evaluate the objective function at sample points and use the provided information to continue 
sampling along promising directions. A prerequisite is an initial point inside the feasible region as well as a set 
of search directions. A number of such methods exist including Nelder-Mead, Hooke-Jeeves, Multi-Coordinate 
Search (MCS); in the current work we use the Hooke-Jeeves method. For completeness we briefly describe the 
method (the reader is referred to [23] for a detailed presentation). 

Hooke- Jeeves is a stencil based method, similar to coordinate descent; however it is a more aggressive search. 
As in all stencil-based methods, a set of search directions v is provided, as well as an array of M scales, s, that 
determine the step length size which the search algorithm is allowed to venture away from the current point when 
investigating for promising search directions. Let us define the design variable vector as x E IR,^, the objective 
function as /(■)> and (without loss of generality) the optimization problem as one of minimizing the objective 
function. In the following pseudocode we present the Hooke-Jeeves method, comprising of the following steps 
[23]: 

• Define: Initial position, xo, search directions, v, and scales, s. 

• search: 

1. Compute f{xo). 
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2. choose scale Si from s. 

3. Exploratory step: 

3.1 Define Xg = xq- 

3.2 Investigate search direction Vj G v: 

— Compute f{xs + SiVj). If f{xs + SiVj) < f{xo) define Xg = Xg + SiVj, move to step 3.3 

— If f{xs + SiVj) > f{xo) compute f{xo - SiVj). 

— If f{xs — SiVj) < f{xo) define Xg = Xg — SiVj, continue to step 3.3. 

3.3 Repeat step 3.2 for j = j + 1. If j = N (all N search directions in v investigated) continue to step 
3.4. 

3.4 Obtain promising direction di = Xg — xq- 

4. If di = 0, then i = i + 1 (reduce scale size to s,+i). If i = M (all M scales in s investigated) terminate 

search iterations. 

5. Pattern move step: 

— Define Xc = xq + 2di. Compute f{xc)- 

- If f{xc) < f{xg), set Xg = Xc- 

6. Set Xq = Xg, and repeat step 1. 

• Optimal point, xq, obtained. 

In Figure 1 we sketch the application of the above pseudocode in a representative optimization problem with 
N = 2; only one iteration is presented. The optimal location in Figure 1 is represented by Xopt, and darker 
contour lines present decreasing values of the objective function. Position xq is chosen to initialize the search 
(step 1), and during the exploratory step (step 3), the objective function is evaluated in the search directions 
Vi and V2 and at a distance si from xq. Two directions are computed to lead to lower objective function values 
(shown with arrows, step 3.2), thus defining Xg at the end of step 3.3 to be at the upper right, and proposing the 
move direction shown with the dotted arrow in Figure 1. During the pattern move step (step 5), an aggressive 
move to Xc is proposed (shown in Figure 1), and f{xc) is computed to be less than f{xg). As a result, the next 
iteration takes place at position Xc, shown in Figure 1 

The method can be easily modified to incorporate inequality constraints for the variables x, by modifying 
the objective function to increase when x moves outside a feasible set. The search algorithm behavior at the 
boundaries is reminiscent of interior point methods. Optimization algorithms that are based only on function 
evaluations are not guaranteed to converge to a global minimum, or even to a minimum. This is due to the fact 
that the necessary optimality conditions in the neighborhood of the result are not computed, and also because a 
search direction may have been neglected. 

Once the optimization algorithm has converged and produced an "apparent best" operating policy profile, 
it may be prudent to restart using the result of the previous search as an initial guess. Restarting the search, 
which invokes again the "large" scales, will cause large perturbations in the operating policy profile, which may 
assist in finding a better local optimum. Once two consecutive runs have produced comparable results, the free 
variable profile is declared "optimal over all scales" [23]. We note that due to the use of KMC simulations, two 
simulation runs will not produce the same value of the objective function for the same policy. This results in 
small variations of the identified policy profiles in two consecutive searches, if small scales are used (scales that 
lead to perturbations in the system state evolution that are comparable to the noise of the KMC simulations). 

4 Numerical Results: NO reduction on Pt surface 

The approach outlined in section 3.1 was initially applied towards the computation of an optimal switching policy 
(between different stationary states) in our simplified NO reduction model. We define a (slightly arbitrary, but 
useful for illustration purposes) objective function in Eq.8 as: 

N 

Q = {k{t) - kgg)^{l - 0.3e-*)T(^,5(f - iT)) 
W = 50[1 - exp{-R{\e{tN) - eggj\ - e))] 
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where e = 0.05, R denotes the ramp function and S denotes the standard Dirac function. We assume that the 
single "manipulated variable" for this problem (the process parameter p of Eq.l) is the (macroscopic) reaction rate 
constant k; this affects the reaction probability in the microscopic simulation. In continuum models the reaction 
rate constant could be manipulated, for example, by changes in temperature, through its Arrhenius temperature 
dependence; our adsorption rate constants could be varied by changes in the gas phase pressure of a species. The 
options that were used for the specific optimization problem arc presented in Table 1. The form of Q reflects the 
formulation of the time-dependent policy through a finite number of decisions, placing a penalty at the time of 
decision. W is a class K function penalizing the deviation of the final state from the desired steady-state heavily 
initially, reaching a plateau at the maximum deviation of 6. 

KMC simulations based on the Gillespie SSA algorithm formed the basis of the coarse time-stepper that was 
used to estimate the coarse-grained system response. A variety of lattice and sample sizes were used in estimating 
the dynamic behavior of the system. Hookc- Jeeves was the algorithm of choice in our search for the optimal 
profile. In Table 2 we present the value of the objective function for the computed optimal parameter profiles, 
through which the effect of the error in the computed optimal profile is implicitly quantified. The lattice size Ni 
was increased with no appreciable change in the expected system response. 

A secondary gain from the lattice size increase was that the ensemble variance d decreased, since d oc (iV;)~^/^. 
Moreover, as the sample size used to compute the coarse response increased in the KMC simulations, the 
ensemble variance d also decreased, as expected, since d oc {Mr)~^^'^. This in turn lead to results from the search 
for the optimal profile of k that are closer to ones obtained when we use the time-stepper of the actual deterministic 
problem (for comparison purposes). The Gillespie stochastic simulation algorithm [15, 17] was chosen so that 
the coarse behavior at large sample sizes can be well approximated through deterministic ODEs, and so that the 
noisy time-stepper optimization results can be compared to the deterministic ones at the appropriate limit. The 
objective value convergence to the computed optimal value from the ODE "direct simulator" comes, of course, 
at the cost of increased computational work. This is only an illustrative example, to validate the approach; it 
would not make sense to use coarse-grained computations when good deterministic equations for the macroscopic 
behavior are known. The equation-free approach is intended for cases where the macroscopic equations are not 
available in closed form. The use of parallel computing can, as we discussed, drastically decrease the necessary 
wall-clock simulation time. In Figure 4 we present the results for Ni = 500 x 500 and Mr = 1000 and compare 
them to the direct macroscopic ODE simulation results. A near-optimal parameter profile is arrived at, due to 
the combination of KMC simulations' noise and the Hookc- Jeeves search algorithm (a search direction that has 
been investigated and characterized as unfavorable is not reinvestigated to conserve CPU time). 

We also used an alternative algorithm (Implicit Filtering, [23]) to compute a near optimal steady state switch- 
ing profile of k{t). The method of implicit filtering uses consecutive boimdcd approximations of the Jacobian 
employing a set (varying at each iteration) perturbation of the design variables and poses a prescribed limit to 
the maximum change of the design variable at each iteration. The time-stepper protocol used had dmax = 0.005, 
with variable lattice size set at Ni = 126 x 126 and Nmax = 4:Ni. Adaptively adjusting Mr, originally set at 
Mr = 252, with Mmax = 4Mr and Mmin = 0.5Mr, the bounds on the ensemble variance were satisfied. We also 
note that a minimum lattice size bound was always enforced in order to ascertain that the discrepancy between 
the expectation of the finite-size KMC simulations and the mean-field ODEs for the same mechanism is negligible 
over the time scales we are working. 

The resulting near optimal time profile of k{t) is presented in Figure 5a, while the near optimal path of NO 
coverage evolution 9{t) is shown in Figure 5b for an averaged realization, and in Figure 5c for a single KMC 
realization with lattice size A*"; = 100 x 100 and reporting horizon 6t = 0.0039. 

In Table 3 we present computational results obtained through incorporating the coarse time-stepper in an 
implicit filtering algorithm. As one might expect, as the variance limit is decreased the algorithm -for a large 
enough lattice- approaches the optimal profile of k{t) computed directly through the coarse ODE model. We also 
observe that a successful scales selection depends heavily on dmax, as scales below a certain limit have an adverse 
effect on the computed near optimal path result (compare the results of the search when the lowest scale is 2~^ 
to the one when it is 2"'^). 

The time-stepper protocol can be combined with the search algorithm to solve the optimization program using 
successively more refined scales. A first search, with large values for the scales, can take place at higher d„iax, 
followed by searches with gradually lower dmax and smaller scales to refine the search for the optimal path; this 
approach may lead to computational savings. When using the KMC with dmax = 0.005 and a lower scale of 2~^ 
for an initial search, followed by a KMC simulation of the system with dmax = 0.0005 and lower scale of 2~^, we 
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find a near optimal path witii v = 10.4110 in a total time of 3981 s. A search using KMC with dmax = 0.001 and 
lower scale of lead to computing a near optimal path of v = 10.4129 in 6858 s (results are shown in Table 3). 

The "optimal" switching path, shown in figure 5b, takes the (expected) phase point from 9ss,s through the 
unstable steady state 9ss,i as shown in Figure 5b. After this is accomplished, we observe that the optimal k{t) 
trajectory rapidly converges back to kss (see Figure 5a). The coarse phase point is now within the region of 
attraction of the steady state O^sj, and no particular switching action is needed to get it there. 

The objective function values presented in this section were computed, for comparison purposes, by integrating 
the coarse system Eq.3 for the optimal switching profile found. During the optimal search using KMC simulations, 
the coarse process model was never used. Several optimization algorithms were explored in conjunction with 
the coarse time-stepper, including Hooke-Jeeves, Nelder-Mead, Implicit-Filtering as well as Multilevel Coordinate 
Search (MCS). Hooke-Jeeves was primarily chosen due to the simplicity of the method and its relative convergence 
speed. Implicit Filtering algorithm was mainly used with the coarse time-stepper to illustrate the variance- 
reduction protocol. Enforcing an upper bound on the variance of the ensemble, along with the selection of the 
value of the scales, provided an order of magnitude estimate of the error in the estimated derivatives during the 
Jacobian approximation; variance reduction techniques may also be useful for this purpose. 

5 Numerical Results: catalytic CO oxidation 

After solving a one-dimensional (coarse grained) example, wc used the same computational approach to find 
optimal switching policies between different stationary states of a simplified CO oxidation model presented in 
section 2.2. The (again slightly arbitrary, illustrative) objective function in Eq.8 was defined as: 

JV 

Q = im - /3,«)2(1 - 0.3e-')T{J2S{t - iT)) 

W = 50[1 - exp{~R{\eco{tN) - eco,.sj\ - e) - i?(|0o.(tjv) - ^o.,.«,/| - e))] 

where e = 0.05, R denotes the ramp function and 5 denotes the Dirac function. The single "manipulated variable" 
for this problem (the process parameter p of Eq.l) was chosen to be the oxygen O2 adsorption rate constant /?, 
which can in principle be manipulated experimentally by varying the oxygen gas phase pressure. The options that 
were used for the specific optimization problem are presented in Table 4. In figure 6 we present the phase portrait 
of the system under the specific nominal process parameter settings. The evolving coarse-grained system has to 
traverse a onc-dimensional separatrix in two phase-space dimensions, in order to enter the basin of attraction of 
the desired CO rich steady-state. 

The kinetic Monte-Carlo time-stepper and the Hooke-Jeeves algorithm were used to compute the coarse 
optimal temporal profile of the O2 adsorption rate (3{t) for large timesteps of T = 0.5 and lattice size of 100 x 100 
and averaging of 200 runs to describe the process evolution. For comparison purposes the optimal profile was 
also identified, using successive quadratic programming (SQP), when the constraints were computed through the 
integration of the mean field deterministic process model of Eq.4 (called the deterministic profile and program, 
respectively, for the rest of the section). The identified temporal profile of /3 is shown in Figure 7a with a blue line 
and is compared to the computed profile of the deterministic program, shown with a green line. The value of the 
objective function is ^'i\/c, 100 x 100.200 = 25.7498, a relative error of 0.55% over the optimal value of vlc = 25.6096 
. The total time for the computation of the optimal profile using KMC (first initial run) was 8 minutes and 16 
seconds. Note that the values of the objective presented -for comparison purposes, and only upon convergence- 
are computed by applying the located optimal policy to the mean field ODE model of Eq.4. Similar results 
were also reached when using the deterministic mean field equations Eq.4 with Nelder-Mead after three search 
initializations and Hooke-Jeeves also after three search initializations. Presenting the CO and O2 coverage in 
Figures 7c and 7d respectively, we observe that the CO oxidation evolves similarly for both (3 profiles and both 
realizations (the deterministic and the KMC-based one). The same observation is made in Figure 7b, the phase 
portrait of the process evolution. We observe that the two evolution profiles, the coarse-grained one using the 
KMC computed optimal policy (green line) and the one computed by integrating Eq.4 lie close to each other at 
all times. 

To obtain a better resolved approximation of the continuous-time optimal manipulation of parameter /3, we 
decreased the time-step to T = 0.1. The discrete-time formulation was subsequently solved using the Hooke-Jeeves 
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algorithm; results are shown in Table 5. Specifically, the search algorithm was reinitialized three times, using the 
previously estimated optimal profile to initialize the search, before an optimal over all scales was declared in the 
case of the noisc-frcc legacy simulator. Similarly, when using the KMC realization of the process, the search was 
terminated when the resulting optimal profile was a slight perturbation to the previous result and the average 
value of the objective (over ten simulations) was within error bounds of the previous value. 

Initially, to solve the optimization problem of Eq.8 with the Hookc- Jeeves search algorithm wc used a KMC 
realization of the process with a Ni = 100 x 100 and Mr = 200 for the coarse time-stepper output; this combination 
identified the /3 temporal profile presented in Figure 8a (green line). We observe that the profile is a perturbation 
around to the resulting f3 profile (blue line) of the deterministic SQP search. 

In Figure 8b we present the phase portrait of the process evolution for the coarse KMC-based HJ computed 
profile (blue line); the system response is closely related to the system behavior under the deterministic, SQP- 
computed, profile (green line). The same observation is also made based on the evolution of CO and O2 coverage, 
shown in Figures 8c and 8d respectively. The closeness of the two phase portraits under the two different switching 
policies suggests that the objective function is relatively insensitive to the change of /3, a further difficulty towards 
the computation of a truly optimal policy. 

To better illustrate the effect of noise originating from the stochastic nature on the time stepper, we run ten 
independent KMC simulations with TV; = 100 x 100, Mr = 200 for the switching policy of Figure 8. The average 
value of the objective function for these ten runs was Vav.MC = 25.1272 with standard deviation (Jmc = 0.0241, 
while the value given by the mean-field ODE model is given at Table 5, 1st line. Based on the intuition gained from 
these runs, we subsequently repeated the third pass of the specific search ten times. The objective function was 
then independently computed for the previously identified optimal profiles using the ODE model (for comparison 
purposes); the average value was computed to be Vav,LC = 25.1928 with standard deviation a^c = 0.0681. 
Let us reiterate that the KMC coarse time-stepper and the ODE time-stepper represent different approximate 
realizations of the same process and as such the predictions of the two time-steppers may lie very close, but in 
general they are not necessarily the same. 

The increase in accuracy comes at a high computational cost, as can be seen from the total CPU time needed 
for each search, when using a KMC realization of the process with a 100 x 100 lattice size and an average of 200 
runs, compared to when using a KMC realization of the process with a 300 x 300 lattice size and an average of 
600 runs. In order to reduce the total computational cost, during the last run (line 3 in Table 5), we initially 
started a search with the accuracy of the coarse timesteppcr implicitly set at a low value (KMC realization with 
100 X 100 lattice size and an average of 200 runs). Once the search converged, we indirectly increased the accuracy 
of the timestepper (KMC realization with 300 x 300 lattice size and an average of 600 runs) and initiated a second 
pass. The computational savings amoimtcd to approximately 55000 seconds. Computational costs can be further 
reduced by adjusting the scales during each new search, taking into account the variance of the coarse time-stepper 
results ensemble. Using a KMC realization of the process with a 300 x 300 lattice size and an average of 600 runs 
for the coarse time-stepper output, the Hooke- Jeeves search algorithm identified the /? temporal profile presented 
in Figure 9a (green line). This compares well with the /3 profile (blue line) resulting from a deterministic SQP 
search. 

Furthermore, in Figure 9b wc present the phase portrait of the process evolution for the coarsc-KMC based, HJ 
computed profile (blue line) and the deterministic SQP-computed profile (green line); we observe that the system 
response is practically the same for both computed profiles of /3{t). The same observation is also made based on 
the evolution of CO and O2 coverage, shown in Figures 9c and 9d respectively, for an averaged realization of the 
system, and in Figures 10a and 10b respectively for a single KMC realization, with time-step 6t = 0.0031. 

The located optimal system trajectory is also shown in Figure 11a for T = 0.5 and Figure lib for T = 0.1, 
with a red line; the blue line represents the separatrix of the system (if it exists) for the current f3{t). The optimal 
switching policy takes the (expected) phase point from 9ss,s through the separatrix of the saddle (unstable) steady 
state 9ss,i; note that during traversing the separatrix line of beta(t) is reduced to a value where no unstable 
steady state (and thus separatrix) exists. After this is accomplished, the optimal /3(t) rapidly converges back to 
(3ss (see Figure 5a). The coarse phase point is now within the region of attraction of the steady state 9ss,f, and 
no particular action is needed. 

Another way to accelerate the compiitation of the optimal switching policy, is by solving a sequence of 
optimization programs, reducing the time-step in each iteration from an initial "coarse" to a final "finer" one. 
In effect the suggestion is to use a method inspired from multigrid optimization methods (e.g., [6, 7]). During 
each iteration, we utilize the result of the previous iteration as an initial guess, which we adapt to the larger 
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number of decisions in time that we need to take. Specifically for this work, we used cubic splines to fit the 
optimal parameter profiles obtained from n-th iteration, and resampled in the new temporal domain to obtain 
the initial guess for the refined search of the n + 1th iteration. Subsequently, during the n + 1th iteration we 
use an off-the-shelf search algorithm to obtain optimal process parameter profiles. As is expected in all multigrid 
methods, the accuracy of the off-the-shelf search algorithm used has to be increased to converge closer to the 
optimal solution. This is done in our case by gradually reducing the size of the scales. Once we have reached 
the "finer" time-step and obtained the final solution using this method, one must investigate if the solution is 
"optimal over all scales" . This is accomplished in our case by performing a new search using scales with a wide 
range of values. In Table 6 wo present the computational efforts needed, following the aforementioned procedure, 
using HJ algorithm with a 100 x 100 lattice with 200 run averaging, KMC realization of the process. In this case 
the scales were kept constant. We observe that we obtain a first pass of the HJ algorithm at T = 0.1 in a total 
of 3983 seconds. 

6 Conclusions 

We presented a computational methodology for the location of coarse-grained optimal operating parameter policies 
for chemically reacting systems described by microscopic/stochastic evolution rules. In particular, we approxi- 
mated optimal operating policies switching bistable reacting systems from one stationary state to another; our 
approach is intended for systems for which macroscopic, coarse evolution equations exist but arc not available in 
closed form. The advantage of the proposed method lies in the establishment, through the coarse time-stepper, 
of a computational bridge between atomistic/ stochastic simulators and traditional (in particular, derivative free) 
optimization algorithms. The approach can be directly extended to systems with higher dimensional expected 
behavior (see for example [33]), and possibly, through matrix- free methods, to systems with infinite dimen- 
sional (spatially distributed, but dissipative) expected behavior [14] . Current efforts are focused on applying this 
methodology to the study of coarse-grained optimal paths associated with transitions and rare events in compu- 
tational chemistry; in this case the coarse timestepper will be estimated from processing the results of an "inner" 
molecular dynamics or an "inner" Monte Carlo simulator ([21, 28]). In this case, the knowledge of appropriate 
macroscopic observables (reaction coordinates) is crucial, and data processing techniques capable of extracting 
such reduced data descriptions (e.g., [42]) become important. The approach we described in this paper can be 
easily combined with "empirical" observables suggested by such data processing algorithms. 
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Table 1: NO reduction process parameters 





Vnlnc 


St(>(Ul\- statcK 


T 


0.25 


ess,s 0.3301 


N 


20 


ess,i 0.6803 


tf 


5 


e^sj 0.9896 




0.45 




a 


1.0 






0.01 
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Table 2: Hooke-Jeeves search results 



Model 


A7 


M, 


()l)jocli\(^ 


r [s] 


Legacy 






10.3709 


129 


KMC 


100 X 100 


200 


10.5418 


674 


KMC 


200 X 200 


400 


10.4155 


4673 


KMC 


300 X 300 


600 


10.3973 


16300 


KMC 


400 X 400 


800 


10.3815 


38469 


KMC 


■500 X 500 


1000 


10.3811 


75307 



* single CPU pentium IV at 2.4 GHz 
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Table 3: Implicit Filtering search results 



AI(k1c1 






Sc< 


lies 




Object 


i ' [s] 


Legacy 












10.3709 


168 


KMC 


0.005 


2" 





•,2 


~3 


10.5461 


297 


KMC 


0.005 


2- 


-0 

1 • 


.,2 


-4 


10.7922 


282 


KMC 


0.0005 


2- 


3 

1 • 


•,2 


-5 


10.4110 


3684 


KMC 


0.001 


2- 


-0 

) • 


•,2 


-4 


10.4129 


6858 



* single CPU pentium IV at 2.0 GHz 
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Table 4: CO oxidation process parameters 





Vnlnc 


Sl('(Ul\- states 


T 


0.25 


(^CO,ss,s 


.13944 


N 


20 


OcO,ss,i 


.67526 




5 


OcO,ss,f 


.97101 


k 


1.0 




.63553 


a 


1.6 


^02,ss,i 


.11452 


7 


0.04 


&02,SS,f 


.00137 


Pes 


3.5 
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Table 5: Hooke Jeeves search results. T = 0.1, t/ = 5 



Model 


Ni 


Mr 


Passes 


Objective 


Total t [s] 


Legacy 






3 


24.8958 


3381 


KMC 


100 X 100 


200 runs 


3 


25.1487 


8182 


KMC 


300 X 300 


600 runs 


2 


24.9997 


57482* 



*First search used a 100 x 100 lattice size KMC realization. 
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Table 6: Hooke Jeeves search results, tf = 5, KMC with 100 x 100 lattice, 200 runs average 



Model 


Timestep 


Objective 


Total t [s] 


KMC 


0.50 


25.7498 


481 


KMC 


0.25 


25.4800 


899 


KMC 


0.10 


25.2734 


2603 


KMC 


0.10* 


25.2607 


3100 



* Second pass at desired timestep. 
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Figure 1: Two-dimensional optimization problem, presenting the steps of one iteration during an optimal search 
using Hooke-Jeeves algorithm. Darker contour lines denote lower values of objective function. 
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Figure 2: Bifurcation diagram of 9 at steady state with respect to k. 
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Figure 3: Bifurcation diagram of 9a and 9b at steady state with respect to (3. 
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Figure 4: Results using Hooke- Jeeves algorithm through numerical integration of Eq.l (blue line) and using KMC 
simulation (green line), a) Optimal temporal profile of process reaction rate k, b) Evolution of NO coverage 9. 
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Figure 5: Results using Implicit Filtering algorithm through numerical integration of Eq.l (blue line) and using 
KMC simulation (green line), a) Optimal temporal profile of process reaction rate k, b) Evolution of NO coverage 
e, c) Evolution of NO coverage 9 for a single KMC reaHzation, Ni = 100 x 100, St = 0.0039. 
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Figure 6: phase portrait of system of Eq.4 presenting the three steady states and the separatrix for /3 
(process parameter values shown in Table 4). 
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Figure 7: Results of Hooke- Jeeves algorithm through numerical integration of Eq.4 (blue line) and using KMC 
simulations with 100 x 100 lattice size and 200 repetitions (green line), a) Optimal temporal profile of O2 adsorption 
rate /3, b) phase portrait of process evolution, c) evolution of CO coverage Oa, d) evolution of O2 coverage Ob 
{tf = 5, TV = 10, I3ss = 3.5). 
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Figure 8: Results of Hooke- Jeeves algorithm through numerical integration of Eq.l (blue line) and using KMC 
simulations with 100 x 100 lattice size and 200 repetitions (green line), a) Optimal temporal profile of O2 adsorption 
rate /3, b) phase portrait of process evolution, c) evolution of CO coverage Oa, d) evolution of O2 coverage Ob 
{tf = 5, TV = 50, Pss = 3.5). 
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Figure 9: Results of numerical integration of Eq.4 (blue line) and KMC simulations with 300 x 300 lattice size 
averaged over 600 runs (green line), a) Optimal temporal profile of O2 adsorption rate /3, identified by Hooke- 
Jeeves search algorithm, b) phase portrait of process evolution, c) evolution of CO coverage 9a, d) evolution of 
O2 coverage 9b [tf = 5, N = 50, f3ss — 3.5). 
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a) ». b) 

Figure 11: Optimal trajectory for system of Eq.4 (red line), identified by Hook- Jeeves search algorithm, and 
evolution of separatrix (blue line), a) for N — 10, b) FOR N — 50, identified hy {tj — 5, Pss — 3.5). 
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